/*
 * Lag.cpp
 *
 *  Created on: Oct 25, 2011
 *      Author: correa
 */

#include <LagSolvers.h>
#include <assert.h>

// Lagrangian relaxation

RepLag::RepLag(Graph * g, Weight<double> * w) : RepStab(g, w) {
}

template <>
RepLag::RepLag(Graph * g, ModelerFactory<RepModeler,Graph> * fm, Weight<double> * w) : RepStab(g, fm, w) {
	subsolver = new CPXLPptr[rm->nreps()];
}

template <>
RepLag::RepLag(Graph * g, ModelerFactory<RepModeler,Graph> * fm, Weight<double> * w, int col, const char * prefix) : RepStab(g, fm, w, col, prefix) {
	subsolver = new CPXLPptr[rm->nreps()];
}

template <>
RepLag::RepLag(Graph * g, ModelerFactory<ARepModeler,Graph> * fm, Weight<double> * w) : RepStab(g, fm, w) {
	subsolver = new CPXLPptr[rm->nreps()];
}

template <>
RepLag::RepLag(Graph * g, ModelerFactory<ARepModeler,Graph> * fm, Weight<double> * w, int col, const char * prefix) : RepStab(g, fm, w, col, prefix) {
	subsolver = new CPXLPptr[rm->nreps()];
}

RepLag::~RepLag() {
	delete subsolver;
}

void RepLag::build(int k) {
	int i;

	nactive = 0;
	const int * cliquerank = rm->cliquerank();
	const int * rd;
	int nreps = rm->nreps();
	for (i = nreps - 1; i >= 0; i--) {
		repm[i] = rm->submodeler(i, 0);
		int nr = repm[i]->graph()->nverts();

		if (nr) {
			int j, cr[nr];
			rd = rm->rd(i, &nr);
			for (j = 0; j < nr; j++)
				cr[j] = cliquerank[rd[j]];
			subsolver[i] = opensolver();
			Stab::build(repm[i], &erowbeg, &erowend, &crowbeg, &crowend, &ub[i], subsolver[i]);

			if (k > 1 || ub[i] > lb) {
				repm[i]->scalerows(erowbeg, erowend-1, 1.0);
				CPXchgrhs (cpxenv, subsolver[i], repm[i]->rcnt(), repm[i]->rmatbeg(), repm[i]->rhs());
				repm[i]->scalerows(crowbeg, crowend-1, 1.0);
				CPXchgrhs (cpxenv, subsolver[i], repm[i]->rcnt(), repm[i]->rmatbeg(), repm[i]->rhs());

				int stat = CPXMIP_FEASIBLE;
				if (k == 1) {
					rm->repedgesymcuts(i);
//					CPXaddrows(cpxenv, subsolver[i], 0, rm->rcnt(), rm->nzcnt(), rm->rhs(), rm->sense(), rm->rmatbeg(), rm->rmatind(), rm->rmatval(), NULL, rm->rowname());
					CPXaddusercuts(cpxenv, subsolver[i], rm->rcnt(), rm->nzcnt(), rm->rhs(), rm->sense(), rm->rmatbeg(), rm->rmatind(), rm->rmatval(), rm->rowname());
					CPXfeasopt(cpxenv, subsolver[i], NULL, NULL, NULL, NULL);
					stat = CPXgetstat(cpxenv, subsolver[i]);
				}

				if (stat == CPXMIP_FEASIBLE) {
					active[nactive++] = i;
					CPXdelmipstarts (cpxenv, subsolver[i], 0, 0);
				}
			}
		}
	}
}

// Max Rep Lags

void MaxRepLag::init(int k) {
	build(k);
	objective();
	int i;
	for (i = 0; i < rm->nreps(); i++)
		repm[i]->delall();
	rm->delall();
}

template <>
MaxRepLag::MaxRepLag(Graph * g, ModelerFactory<RepModeler,Graph> * fm, int k) : RepLag(g, fm, new Weight<double>()), unit(true), hassol(false) {
	init(k);
}

template <>
MaxRepLag::MaxRepLag(Graph * g, ModelerFactory<RepModeler,Graph> * fm, Weight<double> * w, int k) : RepLag(g, fm, w), unit(false), hassol(false) {
	init(k);
}

template <>
MaxRepLag::MaxRepLag(Graph * g, ModelerFactory<RepModeler,Graph> * fm, int col, const char * prefix, int k) : RepLag(g, fm, new Weight<double>(), col, prefix), unit(true), hassol(false) {
	init(k);
}

template <>
MaxRepLag::MaxRepLag(Graph * g, ModelerFactory<RepModeler,Graph> * fm, Weight<double> * w, int col, const char * prefix, int k) : RepLag(g, fm, w, col, prefix), unit(false), hassol(false) {
	init(k);
}

template <>
MaxRepLag::MaxRepLag(Graph * g, ModelerFactory<ARepModeler,Graph> * fm, int k) : RepLag(g, fm, new Weight<double>()), unit(true), hassol(false) {
	init(k);
}

template <>
MaxRepLag::MaxRepLag(Graph * g, ModelerFactory<ARepModeler,Graph> * fm, Weight<double> * w, int k) : RepLag(g, fm, w), unit(false), hassol(false) {
	init(k);
}

template <>
MaxRepLag::MaxRepLag(Graph * g, ModelerFactory<ARepModeler,Graph> * fm, int col, const char * prefix, int k) : RepLag(g, fm, new Weight<double>(), col, prefix), unit(true), hassol(false) {
	init(k);
}

template <>
MaxRepLag::MaxRepLag(Graph * g, ModelerFactory<ARepModeler,Graph> * fm, Weight<double> * w, int col, const char * prefix, int k) : RepLag(g, fm, w, col, prefix), unit(false), hassol(false) {
	init(k);
}

MaxRepLag::~MaxRepLag() {
	if (unit)
		delete w;
}

int MaxRepLag::objective() {
	int i;
	for (i = nactive - 1; i >= 0; i--) {
		int ai = active[i];
		repm[ai]->objective(w);
		CPXchgprobtype (cpxenv, subsolver[ai], CPXPROB_MILP);
		CPXchgobjsen(cpxenv, subsolver[ai], CPX_MAX);
		CPXchgobj(cpxenv, subsolver[ai], repm[ai]->ccnt(), repm[ai]->rmatind(), repm[ai]->rmatval());
		CPXchgctype (cpxenv, subsolver[ai], repm[ai]->ccnt(), repm[ai]->rmatind(), repm[ai]->type());
		hassol = false;
	}
	return 0;
}

int MaxRepLag::solveByEnumeration() {
	hassol = true;
	bestobj = 0.0;

	int i = rm->nreps() - 1;
	int j, nrp;
	while (i >= 0) {
		const int * rp = rm->rp(i, &nrp);
		for (j = 0; j < nrp && S[rp[j]] > 0; j++);
		if (j < nrp)
			i--;
		else
			break;
	}

	if (i >= 0) {
		int si = 0;
		int nrd;
		const int * rd = rm->rd(i, &nrd);
		int SS[nrd];
		for (j = 0; j < nrd; j++)
			SS[j] = S[rd[j]];

		repm[i]->startingpoint(SS);
		rm->rp(i, &nrp);
		printf("starting point has size %d out of %d expected\n",repm[i]->ccnt()+nrp, lb);
		CPXaddmipstarts(cpxenv, subsolver[i], 1, repm[i]->ccnt(), &si, repm[i]->rmatind(), repm[i]->rmatval(), NULL, NULL);
	}

	for (i = 0; i < nactive; i++) {
		solver = subsolver[active[i]];
		rm->rp(active[i], &nrp);
		CPXsetdblparam(cpxenv, CPX_PARAM_CUTLO, bestobj-nrp+1);
		RepLag::solveByEnumeration();
//		printf("stat=%d\n",CPXgetstat(cpxenv, solver));
		if (CPXgetstat(cpxenv, solver) == CPXMIP_OPTIMAL) {
			double repobj;
			CPXgetobjval (cpxenv, solver, &repobj);
			if (repobj + nrp > bestobj)
				bestobj = repobj + nrp;
		}
	}
	return 0;
}

int MaxRepLag::getBest(double * val) {
	*val = bestobj;
	return hassol ? 0 : 1;
}

int MaxRepLag::getNodes(int * nnd) {
	if (hassol) {
		*nnd = 0;
		int i;
		for (i = 0; i < nactive; i++)
			*nnd += CPXgetnodecnt (cpxenv, subsolver[active[i]]);
		return 0;
	}
	return 1;
}
